
# Install and load required packages
# install.packages("bkmr")
library(bkmr)
library(ggplot2)
# 有点看不懂,需要研究下
# # Generate sample data
# set.seed(123)
# n <- 100  # Number of observations
# p <- 5    # Number of predictors
# X <- matrix(rnorm(n * p), nrow = n)  # Predictor matrix
# Y <- rnorm(n) + X[, 1] + X[, 2]^2  # Outcome variable
# 
# # Define kernel specification
# kernel
# k_spec <- kernel.spec("rbfdot", list(sigma = 1))  # Radial basis function kernel
# 
# # Perform BKMR
# bkm_fit <- bkmr(Y = Y, X = X, k_spec = k_spec)
# 
# # Summarize the results
# summary(bkm_fit)
# https://jenfb.github.io/bkmr/overview.html
#https://cran.r-project.org/web/packages/bkmrhat/vignettes/bkmrhat-vignette.html
set.seed(111)
dat <- SimData(n = 50, M = 4)
y <- dat$y
Z <- dat$Z
X <- dat$X
##############################
z1 <- seq(min(dat$Z[, 1]), max(dat$Z[, 1]), length = 20)
z2 <- seq(min(dat$Z[, 2]), max(dat$Z[, 2]), length = 20)
hgrid.true <- outer(z1, z2, function(x,y) apply(cbind(x,y), 1, dat$HFun))

res <- persp(z1, z2, hgrid.true, theta = 30, phi = 20, expand = 0.5, 
             col = "lightblue", xlab = "", ylab = "", zlab = "")

###################
set.seed(111)
fitkm <- kmbayes(y = y, Z = Z, X = X, iter = 10000, verbose = FALSE, varsel = TRUE)
pred.resp.univar <- PredictorResponseUnivar(fit = fitkm)
pred.resp.univar
ggplot(pred.resp.univar, aes(z, est, ymin = est - 1.96*se, ymax = est + 1.96*se)) + 
  geom_smooth(stat = "identity") + 
  facet_wrap(~ variable) +
  ylab("h(z)")
######################################################

pred.resp.bivar <- PredictorResponseBivar(fit = fitkm, min.plot.dist = 1)

ggplot(pred.resp.bivar, aes(z1, z2, fill = est)) + 
  geom_raster() + 
  facet_grid(variable2 ~ variable1) +
  scale_fill_gradientn(colours=c("#0000FFFF","#FFFFFFFF","#FF0000FF")) +
  xlab("expos1") +
  ylab("expos2") +
  ggtitle("h(expos1, expos2)")
##################################################################





set.seed(111)
fitkm <- kmbayes(y = y, Z = Z, X = X, iter = 1000, verbose = FALSE, varsel = TRUE)

ExtractPIPs(fitkm)

TracePlot(fit = fitkm, par = "beta")
TracePlot(fit = fitkm, par = "sigsq.eps")
TracePlot(fit = fitkm, par = "r", comp = 1)

risks.overall <- OverallRiskSummaries(fit = fitkm, y = y, Z = Z, X = X, 
                                      qs = seq(0.25, 0.75, by = 0.05), 
                                      q.fixed = 0.5, method = "exact")

ggplot(risks.overall, aes(quantile, est, ymin = est - 1.96*sd, ymax = est + 1.96*sd)) + geom_pointrange()


